Chapter 8 HMSC analysis

8.1 Load data

load("data/data.Rdata")
load("hmsc/fit_model1_250_10.Rdata")

8.2 Variance partitioning

# Compute variance partitioning
varpart=computeVariancePartitioning(m)

varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
   group_by(variable) %>%
   summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
   tt()

variable	mean	sd
Random: location	37.900015	25.317903
logseqdepth	56.110626	25.796874
sex	4.937460	5.612719
origin	1.051899	1.282563

# Basal tree
varpart_tree <- genome_tree

#Varpart table
varpart_table <- varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(genome=factor(genome, levels=rev(varpart_tree$tip.label))) %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location"))))

#Phyla
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
varpart_tree <- varpart_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

# Add phylum colors next to the tree tips
varpart_tree <- gheatmap(varpart_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
   scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
varpart_tree <- varpart_tree + new_scale_fill()

# Add variance stacked barplot
vertical_tree <-  varpart_tree +
       scale_fill_manual(values=c("#506a96","#cccccc","#be3e2b","#f6de6c"))+
        geom_fruit(
             data=varpart_table,
             geom=geom_bar,
             mapping = aes(x=value, y=genome, fill=variable, group=variable),
             pwidth = 2,
             offset = 0.05,
             width= 1,
             orientation="y",
             stat="identity")+
      labs(fill="Variable")

vertical_tree

8.3 Posterior estimates

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree

# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    #select(genome,sp_vulgaris,area_semi,area_urban,sp_vulgarisxarea_semi,sp_vulgarisxarea_urban,season_spring,season_winter,sp_vulgarisxseason_spring,sp_vulgarisxseason_winter) %>%
    column_to_rownames(var="genome")

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
postestimates_tree <- postestimates_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
      scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add posterior significant heatmap

postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.5, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
        scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
        labs(fill="Trend")

postestimates_tree +
        vexpand(.25, 1) # expand top 

8.3.1 Origin

8.3.1.1 Positively associated genomes with domestic cats

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Positive") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum, class, family, species)%>%
  paged_table()

8.3.1.2 Positively associated genomes feral cats

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Negative") %>%
    arrange(value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum,class,family,species)%>%
  paged_table()

8.3.1.3 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    dplyr::select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    dplyr::select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
  filter(genome %in% estimate$genome) %>%
  arrange(match(genome, estimate$genome)) %>%
  dplyr::select(phylum, colors) %>%
  unique() %>%
  arrange(phylum) %>%
  dplyr::select(colors) %>%
  pull()

inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

8.3.1.4 Correlations

#Compute the residual correlation matrix
OmegaCor = computeAssociations(m)

# Refernece tree (for sorting genomes)
genome_tree_subset <- genome_tree %>%
        keep.tip(., tip=m$spNames) 


#Co-occurrence matrix at the animal level
supportLevel = 0.95
toPlot = ((OmegaCor[[1]]$support>supportLevel)
          + (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean

matrix <- toPlot %>% 
      as.data.frame() %>%
      rownames_to_column(var="genome1") %>%
      pivot_longer(!genome1, names_to = "genome2", values_to = "cor") %>%
      mutate(genome1= factor(genome1, levels=genome_tree_subset$tip.label)) %>%
      mutate(genome2= factor(genome2, levels=genome_tree_subset$tip.label)) %>%
      ggplot(aes(x = genome1, y = genome2, fill = cor)) +
            geom_tile() + 
            scale_fill_gradient2(low = "#be3e2b",
                       mid = "#f4f4f4",
                       high = "#b2b530")+
            theme_void()

htree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

vtree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#create composite figure
grid.arrange(grobs = list(matrix,vtree),
             layout_matrix = rbind(c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1)))

8.3.1.5 Predict responses (origin)

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = c("domestic","feral")
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "origin", 
                      non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
            predict(m, Gradient = ., expected = TRUE) %>%
            do.call(rbind,.) %>%
            as.data.frame() %>%
            mutate(origin=rep(gradient,1000)) %>%
            pivot_longer(!origin,names_to = "genome", values_to = "value")

# weights:  9 (4 variable)
initial  value 101.072331 
final  value 91.392443 
converged

8.3.1.6 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

8.3.1.6.1 Positive associated functions at element level

# Positively associated

unique_funct_db<- GIFT_db[c(2,4,5,6)] %>% 
  distinct(Code_element, .keep_all = TRUE)


element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  paged_table()

#Get community-weighed average GIFTs per sample
genome_counts_row <- genome_counts %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
  column_to_rownames(., "genome") 
#genome_counts_row <- rownames_to_column(genome_counts_row, "genome")
GIFTs_elements_community <- to.community(elements_table,genome_counts_row,GIFT_db)

element_gift_t <- GIFTs_elements_community  %>% 
  t() %>%
#  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

significant_elements<- element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9)

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(4,9)], by = join_by(sample == sample))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift_filt %>%
  select(-sample)%>%
  group_by(Origin)  %>%
  summarise(across(everything(), mean))%>%
   t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

  Elements        Feral         Tame                                         Function
1    B0103 0.7566040000 0.7195779000                Nucleic acid biosynthesis_UDP/UTP
2    D0205 0.1081381000 0.1066891000                Polysaccharide degradation_Pectin
3    D0208 0.2168561000 0.2087181000 Polysaccharide degradation_Mixed-Linkage glucans
4    D0504 0.0096225710 0.0100233660                Amino acid degradation_Methionine
5    D0507 0.1643304000 0.1693777000                   Amino acid degradation_Leucine
6    D0906 0.0001762601 0.0002054726                Antibiotic degradation_Fosfomycin

8.3.1.6.2 Negatively associated functions at element level

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  paged_table()

element_gift_t <- GIFTs_elements_community  %>% 
  t() %>%
#  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

significant_elements<- element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) 

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(4,9)], by = join_by(sample == sample))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift_filt %>%
  select(-sample)%>%
  group_by(Origin)  %>%
  summarise(across(everything(), mean))%>%
   t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

   Elements       Feral        Tame                                         Function
1     B0204 0.546283700 0.595230700                   Amino acid biosynthesis_Serine
2     B0214 0.211179500 0.236935900                Amino acid biosynthesis_Glutamate
3     B0219 0.011526030 0.010872120                     Amino acid biosynthesis_GABA
4     B0302 0.010085548 0.009823828       Amino acid derivative biosynthesis_Betaine
5     B0303 0.233543100 0.236799000       Amino acid derivative biosynthesis_Ectoine
6     B0309 0.076634680 0.071995880    Amino acid derivative biosynthesis_Putrescine
7     B0310 0.055338700 0.046228620    Amino acid derivative biosynthesis_Tryptamine
8     B0401 0.755886800 0.765103900                        SCFA biosynthesis_Acetate
9     B0601 0.333605100 0.356787900             Organic anion biosynthesis_Succinate
10    B0603 0.369766100 0.412079900               Organic anion biosynthesis_Citrate
11    B0708 0.440560800 0.464048500             Vitamin biosynthesis_Cobalamin (B12)
12    B0709 0.004493435 0.004851109 Vitamin biosynthesis_Tocopherol/tocotorienol (E)
13    B0804 0.676473200 0.683577900      Aromatic compound biosynthesis_Dipicolinate
14    D0508 0.038700810 0.048971490                    Amino acid degradation_Lysine
15    D0512 0.198698300 0.225197200                 Amino acid degradation_Histidine
16    D0517 0.039373280 0.046676760                 Amino acid degradation_Ornithine
17    D0601 0.063668170 0.064472310            Nitrogen compound degradation_Nitrate
18    D0603 0.009203820 0.010291730              Nitrogen compound degradation_Urate
19    D0606 0.051570030 0.037713990          Nitrogen compound degradation_Allantoin
20    D0610 0.031016560 0.036298430        Nitrogen compound degradation_Methylamine
21    D0611 0.009801617 0.009503118   Nitrogen compound degradation_Phenylethylamine
22    D0612 0.017850890 0.018892270        Nitrogen compound degradation_Hypotaurine
23    D0801 0.002099865 0.001182548                   Xenobiotic degradation_Toluene
24    D0802 0.002099865 0.001182548                    Xenobiotic degradation_Xylene
25    D0807 0.061055970 0.071058870                  Xenobiotic degradation_Catechol
26    D0816 0.091350430 0.092723640             Xenobiotic degradation_Phenylacetate
27    D0817 0.043241270 0.050879070           Xenobiotic degradation_Trans-cinnamate
28    D0903 0.009801617 0.009503118             Antibiotic degradation_Cephalosporin
29    D0908 0.117656700 0.137340600                 Antibiotic degradation_Macrolide

positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

table <- bind_rows(positive,negative) %>%
  left_join(unique_funct_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) 

table %>%
  mutate(Element=factor(Element,levels=c(table$Element))) %>%
  ggplot(aes(x=mean, y=fct_rev(Element), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.3.1.7 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D02	0.0081357868	-0.003637596	0.020315293	0.825	0.175
B08	0.0074769039	-0.003441048	0.018287846	0.797	0.203
B01	0.0011491533	-0.005813419	0.007937841	0.616	0.384
S01	0.0008741013	-0.010932437	0.012737002	0.563	0.437

# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D08	-1.169835e-03	-0.0022909469	-2.422247e-04	0.035	0.965
B03	-1.067982e-02	-0.0178277310	-3.557808e-03	0.045	0.955
D06	-3.274391e-03	-0.0069404022	-2.258647e-05	0.098	0.902
B04	-8.312913e-03	-0.0195719001	7.637137e-04	0.128	0.872
B06	-7.014874e-03	-0.0169222091	2.281875e-03	0.169	0.831
D07	-1.194568e-02	-0.0300927877	3.894367e-03	0.175	0.825
D03	-4.498865e-03	-0.0130295373	2.733128e-03	0.211	0.789
D05	-1.795499e-03	-0.0075067047	3.157093e-03	0.223	0.777
S03	-9.537274e-03	-0.0320097826	1.427255e-02	0.247	0.753
B02	-3.855277e-03	-0.0125477755	4.170471e-03	0.251	0.749
S02	-4.673357e-03	-0.0146030922	2.873079e-03	0.304	0.696
B07	-3.854268e-03	-0.0156450544	8.129977e-03	0.310	0.690
D09	-1.950214e-03	-0.0080302387	4.401978e-03	0.320	0.680
B09	-3.898566e-06	-0.0005639176	5.023253e-04	0.366	0.634
D01	-3.475224e-04	-0.0051398870	4.039820e-03	0.402	0.598
B10	-2.188595e-05	-0.0003343497	2.311124e-04	0.470	0.530

positive <- function_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- function_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.02,0.02)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.3.2 Sex

8.3.2.1 Negatively associated genomes with male cats

# Compute variance partitioning
varpart=computeVariancePartitioning(m)

varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
   group_by(variable) %>%
   summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
   tt()

variable	mean	sd
Random: location	37.900015	25.317903
logseqdepth	56.110626	25.796874
sex	4.937460	5.612719
origin	1.051899	1.282563

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree
# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    column_to_rownames(var="genome")

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="sexMale", trend=="Negative") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum, class, family)%>%
    tt()

genome	phylum	class	order	family	species	value
bin_Denmark.4	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella stercoris	0.090
bin_Aruba.14	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella stercoris	0.089
bin_Brazil.91	Bacillota	Bacilli	RF39	UBA660	CAG-988 sp003149915	0.082
bin_Brazil.76	Bacillota	Bacilli	RF39	UBA660	CAG-877 sp900554305	0.017
bin_Malaysia.17	Bacillota	Bacilli	RF39	UBA660	NA	0.016
bin_Brazil.45	Bacillota	Bacilli	RF39	UBA660	CAG-877 sp900554315	0.010
bin_Malaysia.16	Bacillota_A	Clostridia	Oscillospirales	Acutalibacteraceae	NA	0.085
bin_Malaysia.78	Bacillota_A	Clostridia	Oscillospirales	Acutalibacteraceae	Ruminococcus_E bromii_B	0.081
bin_Brazil.69	Bacillota_A	Clostridia	Oscillospirales	Acutalibacteraceae	Clostridium_A leptum	0.020
bin_Brazil.83	Bacillota_A	Clostridia	Lachnospirales	Anaerotignaceae	NA	0.069
bin_Denmark.63	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Anaerostipes hadrus	0.091
bin_Spain.11	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.084
bin_Malaysia.3	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.077
bin_Malaysia.45	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia_A wexlerae	0.077
bin_Brazil.1	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.071
bin_Malaysia.7	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	CAG-317 sp900543415	0.070
bin_Denmark.118	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Dorea_A longicatena	0.066
bin_Brazil.3	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900555395	0.065
bin_Malaysia.52	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Fusicatenibacter saccharivorans	0.065
bin_Brazil.89	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900550235	0.062
bin_Denmark.34	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Fusicatenibacter saccharivorans	0.062
bin_Malaysia.73	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.062
bin_Brazil.166	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.057
bin_Brazil.165	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia_A caecimuris	0.055
bin_Denmark.19	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900556835	0.053
bin_Spain.53	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.053
bin_Brazil.54	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Dorea_B phocaeensis	0.051
bin_Brazil.157	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.046
bin_Brazil.32	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Mediterraneibacter torques	0.042
bin_Brazil.8	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia_A sp900547615	0.042
bin_Spain.6	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Schaedlerella glycyrrhizinilytica	0.039
bin_Denmark.109	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Eisenbergiella sp900550285	0.034
bin_Brazil.113	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Anaerobutyricum sp900754855	0.031
bin_Brazil.17	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.030
bin_Spain.37	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	UMGS1472 sp900552095	0.029
bin_Denmark.52	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Ruminococcus_B gnavus	0.023
bin_Brazil.56	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Ruminococcus_A sp000432335	0.022
bin_Spain.67	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Ruminococcus_B sp900544395	0.021
bin_Brazil.28	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.019
bin_Brazil.93	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.019
bin_Malaysia.110	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Catenibacillus sp900557175	0.018
bin_Aruba.43	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Roseburia sp900548205	0.017
bin_Malaysia.108	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.017
bin_Brazil.116	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.016
bin_Malaysia.86	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.015
bin_Brazil.99	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Lachnospira sp900552795	0.013
bin_CaboVerde.61	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Roseburia sp900548205	0.013
bin_Brazil.97	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	UBA9502 sp900538475	0.010
bin_Malaysia.46	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia sp003287895	0.010
bin_Brazil.50	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Lachnospira sp003451515	0.009
bin_Denmark.66	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia sp900120295	0.009
bin_Spain.60	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.009
bin_Malaysia.125	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia sp900539145	0.008
bin_CaboVerde.34	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.005
bin_Aruba.13	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	CAG-81 sp000435795	0.002
bin_Brazil.63	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.002
bin_Denmark.62	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Clostridium_Q sp000435655	0.002
bin_Brazil.105	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.001
bin_Denmark.20	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Enterocloster sp001517625	0.001
bin_Malaysia.98	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.066
bin_Malaysia.9	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Dysosmobacter welbionis	0.065
bin_Aruba.28	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.044
bin_Brazil.177	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.044
bin_Malaysia.116	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Flavonifractor plautii	0.043
bin_Malaysia.34	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Lawsonibacter sp000177015	0.040
bin_Brazil.137	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.027
bin_Aruba.52	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.094
bin_Aruba.31	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	Negativibacillus sp000435195	0.080
bin_Malaysia.102	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	NA	0.080
bin_Malaysia.30	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	Phocea massiliensis	0.055
bin_Denmark.72	Bacillota_C	Negativicutes	Selenomonadales	Selenomonadaceae	Megamonas funiformis	0.082
bin_Denmark.85	Bacillota_C	Negativicutes	Selenomonadales	Selenomonadaceae	NA	0.077
bin_Brazil.5	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella sp900540415	0.080
bin_Malaysia.18	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotellamassilia sp000437675	0.051
bin_Brazil.48	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides sp900766005	0.050
bin_Malaysia.117	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotellamassilia sp900541335	0.047
bin_CaboVerde.18	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella copri	0.038
bin_Aruba.10	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella sp900544825	0.037
bin_Malaysia.64	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Paraprevotella clara	0.033
bin_Brazil.163	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella lascolaii	0.029
bin_Spain.4	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola sp000436795	0.027
bin_Spain.48	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides fragilis	0.019
bin_Brazil.96	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides uniformis	0.014
bin_Malaysia.105	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola massiliensis	0.014
bin_Spain.21	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola coprophilus	0.012
bin_Malaysia.121	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides stercoris	0.011
bin_Brazil.43	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola vulgatus	0.010
bin_Denmark.30	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola sp900546645	0.010
bin_Brazil.103	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola plebeius_A	0.009
bin_Malaysia.77	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola sp900542985	0.006
bin_Brazil.6	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola coprocola	0.005
bin_Brazil.110	Bacteroidota	Bacteroidia	Bacteroidales	Barnesiellaceae	Barnesiella intestinihominis	0.024
bin_Malaysia.13	Bacteroidota	Bacteroidia	Bacteroidales	Muribaculaceae	NA	0.075
bin_Spain.23	Bacteroidota	Bacteroidia	Bacteroidales	Muribaculaceae	CAG-279 sp900541935	0.051
bin_CaboVerde.37	Bacteroidota	Bacteroidia	Bacteroidales	Porphyromonadaceae	NA	0.058
bin_Brazil.11	Bacteroidota	Bacteroidia	Bacteroidales	Rikenellaceae	Alistipes putredinis	0.045
bin_Malaysia.131	Bacteroidota	Bacteroidia	Bacteroidales	Rikenellaceae	Alistipes putredinis	0.043
bin_Brazil.38	Bacteroidota	Bacteroidia	Bacteroidales	Rikenellaceae	Alistipes communis	0.039
bin_Brazil.86	Bacteroidota	Bacteroidia	Bacteroidales	Rikenellaceae	Alistipes dispar	0.030
bin_Brazil.138	Bacteroidota	Bacteroidia	Bacteroidales	Tannerellaceae	Parabacteroides sp000436495	0.023
bin_Brazil.124	Bacteroidota	Bacteroidia	Bacteroidales	Tannerellaceae	NA	0.019
bin_Brazil.160	Bacteroidota	Bacteroidia	Bacteroidales	Tannerellaceae	Parabacteroides distasonis	0.013
bin_CaboVerde.10	Campylobacterota	Campylobacteria	Campylobacterales	Campylobacteraceae	NA	0.070
bin_Brazil.68	Campylobacterota	Campylobacteria	Campylobacterales	Campylobacteraceae	Campylobacter_D helveticus	0.049
bin_Spain.26	Campylobacterota	Campylobacteria	Campylobacterales	Campylobacteraceae	Campylobacter_D upsaliensis	0.047
bin_Brazil.145	Cyanobacteria	Vampirovibrionia	Gastranaerophilales	Gastranaerophilaceae	NA	0.044
bin_Brazil.140	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_A sp900543175	0.029
bin_Malaysia.6	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_B sp900545035	0.015
bin_Brazil.159	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_B sp900541465	0.012
bin_Brazil.146	Pseudomonadota	Alphaproteobacteria	RF32	CAG-239	CAG-495 sp001917125	0.074
bin_Malaysia.50	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	NA	0.084
bin_Brazil.82	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum succiniciproducens	0.046
bin_CaboVerde.33	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum sp900543125	0.045
bin_CaboVerde.55	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum_A thomasii	0.036
bin_Brazil.111	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Succinivibrio sp000431835	0.024

8.3.2.2 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="sexMale") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    dplyr::select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="sexMale") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    dplyr::select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
  filter(genome %in% estimate$genome) %>%
  arrange(match(genome, estimate$genome)) %>%
  dplyr::select(phylum, colors) %>%
  unique() %>%
  arrange(phylum) %>%
  dplyr::select(colors) %>%
  pull()

inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

8.3.2.3 Predict responses

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = c("Male","Female","Unknown")
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "sex", 
                      non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
            predict(m, Gradient = ., expected = TRUE) %>%
            do.call(rbind,.) %>%
            as.data.frame() %>%
            mutate(sex=rep(gradient,1000)) %>%
            pivot_longer(!sex,names_to = "genome", values_to = "value")

# weights:  4 (3 variable)
initial  value 63.769541 
final  value 61.728471 
converged

8.3.2.4 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- pred %>%
  group_by(sex, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "sex") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="sex")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "sex") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

8.3.2.4.1 Positive associated functions at element level

# Positively associated

unique_funct_db<- GIFT_db[c(2,4,5,6)] %>% 
  distinct(Code_element, .keep_all = TRUE)


element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)

# A tibble: 27 × 9
   trait    mean       p1      p9 positive_support negative_support Domain       Function                           Element       
   <chr>   <dbl>    <dbl>   <dbl>            <dbl>            <dbl> <chr>        <chr>                              <chr>         
 1 B0212 0.0295  0.00137  0.0554             0.909            0.091 Biosynthesis Amino acid biosynthesis            Arginine      
 2 B0220 0.00838 0.000520 0.0159             0.909            0.091 Biosynthesis Amino acid biosynthesis            Beta-alanine  
 3 B0310 0.0189  0.00558  0.0325             0.959            0.041 Biosynthesis Amino acid derivative biosynthesis Tryptamine    
 4 B0303 0.0146  0.000861 0.0264             0.913            0.087 Biosynthesis Amino acid derivative biosynthesis Ectoine       
 5 B0307 0.0408  0.00190  0.0709             0.906            0.094 Biosynthesis Amino acid derivative biosynthesis Spermidine    
 6 B0901 0.00110 0.000165 0.00258            0.917            0.083 Biosynthesis Metallophore biosynthesis          Staphyloferrin
 7 B0104 0.0328  0.00473  0.0561             0.925            0.075 Biosynthesis Nucleic acid biosynthesis          CDP/CTP       
 8 B0105 0.0245  0.00614  0.0419             0.924            0.076 Biosynthesis Nucleic acid biosynthesis          ADP/ATP       
 9 B0106 0.0139  0.000246 0.0265             0.901            0.099 Biosynthesis Nucleic acid biosynthesis          GDP/GTP       
10 B0605 0.0618  0.0328   0.0910             0.966            0.034 Biosynthesis Organic anion biosynthesis         D-lactate     
# ℹ 17 more rows

8.3.2.4.2 Negatively associated functions at element level

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)

# A tibble: 19 × 9
   trait       mean         p1          p9 positive_support negative_support Domain       Function                      Element           
   <chr>      <dbl>      <dbl>       <dbl>            <dbl>            <dbl> <chr>        <chr>                         <chr>             
 1 B0215 -0.0600    -0.0910    -0.0281                0.028            0.972 Biosynthesis Amino acid biosynthesis       Histidine         
 2 B0216 -0.0477    -0.0710    -0.0229                0.031            0.969 Biosynthesis Amino acid biosynthesis       Tryptophan        
 3 B0208 -0.0286    -0.0528    -0.000896              0.093            0.907 Biosynthesis Amino acid biosynthesis       Valine            
 4 B0209 -0.0274    -0.0539    -0.0000683             0.1              0.9   Biosynthesis Amino acid biosynthesis       Isoleucine        
 5 B1012 -0.00822   -0.0129    -0.00371               0.007            0.993 Biosynthesis Antibiotic biosynthesis       Fosfomycin        
 6 B0704 -0.0585    -0.0919    -0.0262                0.019            0.981 Biosynthesis Vitamin biosynthesis          Pantothenate (B5) 
 7 B0711 -0.0331    -0.0528    -0.0143                0.041            0.959 Biosynthesis Vitamin biosynthesis          Menaquinone (K2)  
 8 B0710 -0.0177    -0.0297    -0.00521               0.059            0.941 Biosynthesis Vitamin biosynthesis          Phylloquinone (K1)
 9 B0703 -0.0134    -0.0263    -0.00142               0.077            0.923 Biosynthesis Vitamin biosynthesis          Niacin (B3)       
10 D0516 -0.0000411 -0.0000755 -0.00000151            0.008            0.992 Degradation  Amino acid degradation        Beta-alanine      
11 D0503 -0.000644  -0.00129   -0.0000797             0.056            0.944 Degradation  Amino acid degradation        Cysteine          
12 D0906 -0.00350   -0.00696   -0.000488              0.04             0.96  Degradation  Antibiotic degradation        Fosfomycin        
13 D0613 -0.0604    -0.0970    -0.0233                0.038            0.962 Degradation  Nitrogen compound degradation Taurine           
14 D0209 -0.0461    -0.0729    -0.0201                0.017            0.983 Degradation  Polysaccharide degradation    Xylans            
15 D0205 -0.0200    -0.0312    -0.00907               0.02             0.98  Degradation  Polysaccharide degradation    Pectin            
16 D0210 -0.0160    -0.0316    -0.00271               0.056            0.944 Degradation  Polysaccharide degradation    Beta-mannan       
17 D0202 -0.0170    -0.0300    -0.00368               0.064            0.936 Degradation  Polysaccharide degradation    Xyloglucan        
18 D0201 -0.0194    -0.0323    -0.00546               0.069            0.931 Degradation  Polysaccharide degradation    Cellulose         
19 D0203 -0.0154    -0.0284    -0.00110               0.093            0.907 Degradation  Polysaccharide degradation    Starch            

positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
#      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
#      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.3.2.5 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- pred %>%
  group_by(sex, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "sex") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="sex")
   })

#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "sex") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D07	0.0370879175	1.865579e-02	0.054600433	0.990	0.010
B03	0.0132361030	2.560580e-03	0.022378474	0.924	0.076
B01	0.0112990801	1.536807e-03	0.021070397	0.922	0.078
B09	0.0008303995	7.280411e-06	0.001905481	0.901	0.099
D06	0.0042811880	-2.679345e-05	0.009422817	0.898	0.102
D03	0.0081369808	-1.021315e-03	0.017939832	0.879	0.121
B04	0.0089402186	-1.852843e-03	0.019829339	0.870	0.130
D09	0.0045467180	-4.073236e-03	0.012796230	0.819	0.181
D05	0.0033674482	-7.647938e-03	0.011635638	0.809	0.191
B06	0.0097564043	-5.525707e-03	0.025204295	0.807	0.193
S03	0.0152479121	-1.896561e-02	0.043166598	0.799	0.201
D08	0.0003058347	-3.947465e-04	0.001036747	0.685	0.315
D01	0.0004439636	-5.283954e-03	0.006055686	0.646	0.354

# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
B10	-0.0006188166	-0.0009435699	-0.0002871661	0.029	0.971
D02	-0.0148823215	-0.0276192261	-0.0021080715	0.074	0.926
B07	-0.0093042984	-0.0214646905	0.0021940396	0.134	0.866
S01	-0.0128671477	-0.0284627187	0.0050748248	0.139	0.861
S02	-0.0033662258	-0.0116276565	0.0056703820	0.179	0.821
B02	-0.0035628916	-0.0153096303	0.0077135414	0.365	0.635
B08	-0.0024022736	-0.0160211150	0.0099821510	0.514	0.486

positive <- function_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- function_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
#      xlim(c(-0.02,0.02)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")